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This paper presents an accurate numerical method for solving a class of fractional optimal control problems 
(FOCPs). In the proposed technique, we approximate FOCPs and end up with a finite dimensional problem. The method is 
based upon the useful properties of Chebyshev polynomials approximation and Rayleigh-Ritz method to reduce FOCPs to 
solve a system of algebraic equations. Also, we use an approximated formula of the Caputo fractional derivative. Special 
attention is given to study the convergence analysis and derive an error upper bound of the proposed approximate formula. 
Illustrative examples are included to demonstrate the validity and the applicability of the proposed technique. 
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Fractional derivatives have recently played a significant role in many areas of sciences, engineering, fluid 
mechanics, biology, physics and economies. Some numerical methods for solving fractional differential equations were 
appeared in (Frederico et al. (2008)-Khader et al. (2013), Sweilam et al. (2007)-Sweilam et al. (2014)). Many real-world 
physical systems display fractional order dynamics that is their behavior is governed by fractional order differential 
equations. For example, it has been illustrated that materials with memory and hereditary effects, and dynamical processes, 
including gas diffusion and heat conduction, in fractal porous media can be more adequately modeled by fractional order 
models than integer order models. 

The general definition of an optimal control problem requires the minimization of a criterion function of the states 
and control inputs of the system over a set of admissible control functions. The system is subject to constrained dynamics 
and control variables. Additional constraints such as final time constraints can be considered. In (Agrawal (2004)), 
the author introduced an original formulation and a general numerical scheme for a potentially almost unlimited class of 
FOCPs. FOCP is an optimal control problem in which the criterion and/or the differential equations governing the 
dynamics of the system contain at least one fractional derivative operator. A general formulation and a solution scheme for 
FOCPs were first introduced in (Agrawal (2004)) where fractional derivatives were introduced in the Riemann-Liouville 
sense, and FOCP formulation was expressed using the fractional variational principle and the Lagrange multiplier 
technique. The state and the control variables were given as a linear combination of test functions, and a virtual work type 
approach was used to obtain solutions. 
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Representation of a function in terms of a series expansion using orthogonal polynomials is a fundamental 
concept in approximation theory and form the basis of the solution of differential equations (Khader et al. (2013), Sweilam 
et al. (2010)). Chebyshev polynomials are widely used in numerical computation. One of the advantages of using 
Chebyshev polynomials as a tool for expansion functions is the good representation of smooth functions by finite 

x(t) 

Chebyshev expansion provided that the function v ' is infinitely differentiable. The coefficients in Chebyshev expansion 
approach zero faster than any inverse power in n as n goes to infinity. Khader (2011), introduced a new approximate 
formula of the fractional derivative and used it to solve numerically the fractional diffusion equation. 

In the present paper, we focus on optimal control problems with the quadratic performance index and the dynamic 
system with the Caputo fractional derivative. Our tools for this aim are the shifted Chebyshev ortho normal basis and the 
collocation method. 

We implement the proposed algorithm for solving the following FOCP in the form 

1 1 9 2 

minimum: / = — \[p{t)x (t) + q(t)u (t)]dt, 

20 (1) 
subject to the system dynamics 

D^ a \{t) = a(t)x(t) + b{t)u(t), (2) 

x(0) = x n 0 

subject to the initial condition u . Where u ^ u refers to the order of the Caputo fractional derivatives, 

pit) > 0, q(t) > 0, a(t) ^ 0 and b(t) arg gj ven f unct j ons Many authors studied these problems with different 
numerical methods. For more details about these problems, see (Agrawal (2004) - Arrow et al. (1958), Frederico et al. 
(2008)). 

PRELIMINARIES AND NOTATIONS 

In this section, we present some necessary definitions and mathematical preliminaries of the fractional calculus 
theory required for our subsequent development. 

The Caputo Fractional Derivatives 
Definition 1 

The Caputo fractional derivative operator u of order u is defined in the following form 

D icQ m= _}_) f (m) (*\ dt , a>0 , 
Y{m-a)h(t-T) a - m+X 

where" 1 " 1 <flr - m ' mG N ' t>0 ' 

Similar to integer-order differentiation, Caputo fractional derivative operator is a linear operation 
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D (a) (Ag(t) + juh(t)) = AD (a) g(.t) + iuD {a) h{t), (3) 

where an " are constants. For the Caputo's derivative we have — ' C is a constantand 



D (a) t n = < 



0, for neNg and n<|~flf|, 

nn + D t n-a for n e n and n > \a\ . 
T(n + \-a) 0 1 1 



(4) 



\(X~\ ry 

We use the ceiling function 1 1 to denote the smallest integer greater than or equal to " and 



N n ={0,1,2,...}. aeN 

u Recall that for ' the Caputo differential operator coincides with the usual differential operator 

of integer order. For more details on fractional derivatives definitions and its properties see (Lotfi et al. (2011), Podlubny 
(1999)). 

The Definition and Properties of the Shifted Chebyshev Polynomials 

The well-known Chebyshev polynomials are defined on the interval [-1,1] and can be determined with the aid of 
the following recurrence formula (Snyder (1966), Sweilam et al. (2007)) 

T n+ l(z) = 2zT n (z)-T n _ l (z), T Q (z) = l, T l (z) = z, n = l, 2,.... 

T (z) 

The analytic form of the Chebyshev polynomials n of degree n is given by 

t ( \ L v ( Wnn~2i-\ n(n-i-\)\ n -2i 0 , 

T n (z)= Z (-1) 2 z , n = 2,3,... 

H i=0 (i)\(n-2i)\ (5) 

where [^^] denotes the integer part of n l 2 

In order to use these polynomials on the interval [0,1], we define the so called shifted Chebyshev polynomials by 
introducing the change of variable Z — 2t -\ 
Lemma 1 

* 

T (t) 

The analytic form of the shifted Chebyshev polynomials n of degree n is given by 



" * ^n-k n 2k n(n + k-l)\ l 0 , 

T n (t)= I (-1) 2 — — —t , n = 2,3,. 

k=0 (2k)\(n-k)\ 



(6) 



Proof 



Tn(t) = T 0 (yTt), 



Since the shifted Chebyshev polynomials can be defined by Ln then by substituting in Eq. (5) 
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we can obtain 

T *(t\ - 1„ v ( W^ln-li-l {2n-i-\)\ n -i „_„ 

T n (t) = 2n I '(-1) 2 — — — t , n-2,3,... 

i=0 (i)\(2n-2i)\ 

Now, if we put k — n — i m tne a b ove equation we obtain the desired result (6). 



Note that from Eq. (6), we can see that n 



T n {Q) = (-\) n , and T n {\) = \ 



DERIVATION AN APPROXIMATE FORMULA FOR FRACTIONAL DERIVATIVES USING 
CHEBYSHEV 

Series Expansion 

x(t) 

The function v ' , which belongs to the space of square integrable functions in [0, 1], may be expressed in terms 
of shifted Chebyshev polynomials as 



(7) 



OO ^ 

x(t)= I c-T. (0, 

c- 

where the coefficients ' are given by 

l]x(t)T*(t) 2]x(t)TAt) . 

c n =— — , af, c,-=— — , af, i=l,2,.„. 

In practice, only the first (m+l)-terms of shifted Chebyshev polynomials are considered. Then we have 

x m (t)= I c-7j. (r). 

*=0 (9) 

x (t) 

The main approximate formula of the fractional derivative of m is given in the following theorem. 
Theorem 1 

(Khader (2011)) 

Let X ^ be approximated by Chebyshev polynomials as (9) and also suppose & ^ ^ , then 

1111 (10) 
where ls g,ve„b y (i-t)!(2t)ir(* + l-a) 

Also, in this section, we will state some theorems to study the convergence of the proposed approximate formula 
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and derive an upper bound of the error. 
Theorem 2 

(Doha et al. (2011)) 

The Caputo fractional derivative of order for the shifted Chebyshev polynomials can be expressed in terms of 
the shifted Chebyshev polynomials themselves in the following form 



D (a) (T*(t))= I 



fc=fa| 7 



I e- ik Tj (t), 

,_Q l,J,K 1 



(11) 



0. 

i,J,k 



where 
Theorem 3 

(Sweilam et al. (2013)) 



(-1)' k 2i(i + k-iy.T(k-a+^) 

h T(k + -k)\T(k-j-a + l)T(k + j-a+l) 
J 2 



h Q =l, hj=2, 7 = 2,3,... . 



The error 



Ej(m) 



D^ a) x{t)-D^ a) x m {t) 



. , D^x(t). D (a) x m (t)., , 
in approximating v ' by lsbounded by 



< 



z=m+l 



Z Z 0; 

*=fa] 7=0 



(12) 



PROCEDURE SOLUTION USING CHEBYSHEV COLLOCATION METHOD 



In this section, we implement the proposed algorithm with using the presented approximate formula of fractional 
derivative (10) to solve numerically the proposed problem of FOCPs defined in (l)-(2). 

The procedure of the implementation is given by the following steps 

Substitute by Eq.(2) into Eq.(l) gives (Lotfi et al. (201 1)) 



minimum: / = - }[p(t)x 2 (t) + [D^x(t) - a(t)x(t)] 2 ]dt. 



26 



(b(t)Y 



(13) 



x(t) 



Approximate the function v ' using Chebyshev polynomial expansion in formula (9) andits 

D ((X) x(t) 

Caputo fractional derivative v ' using the proposed approximate formula (10), then 

Eq.(13) transformed to the following approximated form 
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minimum J =-\ [p(t)( I c-T- X S c-w}? J r a - a(f) I c-T. (t)f]dt, (14) 

2 0 i=0 ' 1 (b(t)) 2 i=\a]k^a] 1 lJc i=0 1 ' 

(a) 

where is defined in (10). 

The integral term in Eq.(14) can be found using composite trapezoidal integration technique as 

r m * 9 ait) m i (a) k-rv m * ? 

\[p(t)( Z cT. (O) 2 +-^y[ I I c f w}" >t k a - a(t) Z c.7. (t)] Z ]dt 

I i=0 1 1 {b{t)) 2 i^ a -\k^a\ 1 hk i=0 1 1 

h N-l 
= -[Q(t 0 ) + Q(t N ) + 2 Z Q(t k )l 
I k=l 

m * ? q(t) m i (a) k-a m * 1 

wt)=p(tx z c.t. (or +-^-[ z z c i w it ^ - a ^ s c,t;. (or, 



(15) 



1 

, for an arbitrary integer N, 1 . So, we can write Eq.(15) in thefollowing 

approximated form 



h ~ N ,. . „ t-=ih, I = 0,1,..., AT 



h N-l 
J(c Q ,c v ...,c m ) = -[a(t Q ) + a(t N ) + 2 z ^)]- 
4 & — 1 



(16) 

The extremal values of functional of the general form (16), according to Rayleigh-Ritz method 

gives 

K =0 *=o ±L=o 

so, after using the boundary conditions, we obtain a system of algebraic equations. 

Cn;Cl; '''TTI X(t) 

Solve the algebraic system to obtain u 1 , then the function v ' which extremes FOCPs 

has the following form 

x m (t) = Z c-T* {t), u(t)=-^-(D (a) x(t)-a(t)x(t)). 

(17) 

APPLICATIONS AND NUMERICAL RESULTS 

In this section, we demonstrate the capability of the introduced approach. To achieve thisaim, we solve two 
widely used examples from the literature. The introduced problems are statedin the traditional FOCP framework and then 
reformulated via our introduced methodology. 



Impact Factor (JCC): 4.2949 



Index Copernicus Value (ICV): 3. 0 



An Approximate Solution for Fractional Optimal Control 71 
Problems Using Chebyshev Pseudo-Spectral Method 

Problem 1 (Linear time-invariant problem) 

Consider the following linear time invariant problem, which described by the following fractional optimal control 
problem (Agrawal (2004) - Agrawal et al. (2007)) 

1 1 o 2 
minimum: / = — j [x (?) + u (t)]dt, 

2 0 (18) 
subject to the system dynamics 

(a) 

D K J x{t) = -x(t) + u{t), 0<or<l, (19) 



subject to the initial condition ^ . 



Our aim is to find the control u ^ which minimizes the quadratic performance index ^ . Fortius problem we 
have the exact solution in the case of ^ ~ 1 as follows 

x(t) = cosh(V20 + fi sinh(V2r), u(t) = (1 + yfcp) cosh(V2t) + (Jl + p) sinh( sflt), 

^ = cosh(V2)+V2sinh(V2) ^ Q9g 
where V2cosh(V2)+sinh(V2) ~ 

We will implement the proposed algorithm described in the previous section with m=3 andapproximate the 
x(t) 

solution v ' as follows 

3 * 
x m (t)= I c-T. (t). 
i=0 1 

The procedure of the implementation is given by the following steps 
Substitute by Eq.(19) into Eq.(18) gives 
1 



(20) 



J =-l[x 2 (t) + (D^x(t) + x(t)f]dt. 
2 0 

x(t) 

Approximate the function v ' using Chebyshev polynomials expansion in formula (9) and its 

D (a) x(t) 

Caputo fractional derivative v ' using the proposed approximated formula (10), then 

Eq.(21) transformed to the following form 

1 1 3 * 9 3 i (sy\ U_fy 3 * 9 

J = - [( I c.T. (0) 2 + [Z I a + Z c.T. (tT]dt, 

2q z=0 1 1 i=U=l 1 l > k i=0 1 1 



(21) 



(22) 
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(a) 

Ik 

where ' is defined in (10). 

The integral term in Eq.(22) can be found using composite trapezoidal integration technique as 

13 # 9 3 i (fy\ 1 r _r, 3*9 h N—l 

[( Z c.t. (or +[Z I c-wW a + z c-r. (or]*=-[n(^ n )+a(^)+2 z ft(f, )], (23) 

0 j=0 z 1 i=\k=\ a i=0 11 2 u « fe=1 * 

fl(0 = ( 2 cT. (0) Z +[Z Z W? J r +2 c i T i (')] , 
/=0 z z f=U=l z z ' /c /=0 z z 



»=1 



" , for an arbitrary integer iV ' 1 
approximated form 



,t-=ih, i = 0,l,..,,N 



So, we can write Eq.(22) in thefollowing 



h N-l 

j(c 0 ,c v c 2 ,c 3 )=-[a(t Q )+a(t N )+2 z &(t k )i 

4 & — 1 

The extremal values of functional of the general form (24), according to Rayleigh-Ritz method gives 

3L =0t K =0t *L =0t *L=o, 

3cq 

Solve the algebraic system using the Newton iteration method to obtain , then the 

function v y which extreme FOCPs (19) has the following form 



(24) 



x m (t) = Z c.T* {t\ u(t) = D^ a) x{t) + x(t). 
?'=0 



(25) 



■ 1 1 1 r 



Exact solution atr/=l 

♦ Approximate solution at a=1 

Approximate solution atu=0.9|5 
□ Approximate solution at a=0.8 i 

Approximate solution at«=0.75 




Figure 1: The Behavior of X ^ for Problem 1 at m = 3 with Different Values of & . 
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Approximate solution at m=4 

+ Approximate solution at m=6 




Figure 2: The Behavior of X ^ for Problem 1 at ^ = 0-8 w ith m=4 and m =6 

The behavior of the numerical solutions of this problem with different values of and m are given in figures 1- 

2. Where in figure 1, the numerical solution X ^ at m=3 for different values of & an( j me exact solution at ^ ~ 1 are 

plotted. In figure 2, the numerical solutions a j CC — 0.8 f or different values of m(m = 4,6) are plotted. 

The solution obtained using the suggested method is in excellent agreement with the already exact solution and 
show that this approach can be solved the problem effectively. It is evident that the overall errors can be made smaller by 
adding new terms from the series (20). Comparisons are made between approximate solutions and exact solutions to 
illustrate the validity and the great potential of the proposed technique. 

Also, from our numerical results we can see that these solutions are in more accuracy of those obtained in 
(Agrawal (2008), Lotfi et al. (201 1)). 

Problem 2 (Linear time-variant problem) 

In this example, we consider the linear time variant fractional optimal control problem (Agrawal (2008), Agrawal 
et al. (2007)) 



1 1 7 2 
minimum: / = — J [jc (t) + u (t)]dt, 

20 



(26) 



subject to the system dynamics 



(27) 



D^x(t) = t x(t) + u(t), 0<a<l, 

x(0) = 1 
subject to the initial condition v ' 

Our aim is to find the control U ^ which minimizes the quadratic performance index J. 
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We will implement the proposed algorithm as described in the previous section with 



m = 3 



0.95 




0.75 



0.65 



Figure 3: The Behavior of X ^\or Problem 2 at m ~ 3 with Different Values of a 



= 0.8 




Figure 4: The Behavior of X ^ for Problem 2 at a ~ 0-8 w i th m - 4 and m - 6. 

The behavior of the numerical solutions of this problem with different values of m and ^ are given in figures 3 

and 4. Where in figure 3, the numerical solutions X ^ at m =3 for different values of ^ and the exact solution at ^ ~ ^ 

are plotted. In figure 4, the numerical solutions X ^ at ^ — 0.8 f or different values of m (in = 4, 6) are plotted. 

CONCLUSIONS AND REMARKS 

In this article, we introduced an accurate numerical scheme for solving a wide class of fractional optimal control 
problems. In the proposed method, FOCP is transformed to solve a system of algebraic equations. The error upper bound 
of the proposed approximate formula is stated and proved. The results show that the algorithm converges as the number of 
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m terms is increased. The solution is expressed as a truncated Chebyshev series and so it can be easily evaluated for 
arbitrary values of t using any computer program without any computational effort. From illustrative examples, it can be 
seen that this approach can obtain very accurate and satisfactory results. For all examples, the solution for the integer order 
case of the problem is also obtained for the purpose of comparison. In the end, fromour numerical results using the 
proposed method, we can conclude that, the solutions are in excellent agreement with the exact solution and better than the 
numerical results obtained in Agrawal and Lotfiapproaches. All computational calculations are made by Matlab program 8. 
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